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Abstract 

We discuss lattice simulations of the ground state of dilute neutron matter at next- 
to-leading order in chiral effective field theory. In a previous paper the coefficients of the 
next-to-leading-order lattice action were determined by matching nucleon-nucleon scattering 
data for momenta up to the pion mass. Here the same lattice action is used to simulate the 
ground state of up to 12 neutrons in a periodic cube using Monte Carlo. We explore the 
density range from 2% to 8% of normal nuclear density and analyze the ground state energy 
as an expansion about the unitarity limit with corrections due to finite scattering length, 
effective range, and P-wave interactions. 
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I. INTRODUCTION 



This is the second of a pair of papers studying chiral effective field theory on the lattice 



at next-to- leading order. In the first paper 



l| we used nucleon-nucleon scattering data at 



low energies to determine unknown operator coefficients of the next-to-leading-order lattice 
action. We also tested model independence of the effective theory at fixed lattice spacing by 
computing next-to-leading-order corrections for two different leading-order lattice actions. 
In this paper we use the Gaussian-smeared lattice actions LO2 and NLO2 defined in [l| to 
simulate dilute neutron matter in a periodic cube. We probe the density range from 2% to 
8% of normal nuclear matter density. Neutron-rich matter at this density is likely present 
in the inner crust of neutron stars 0, sl. The Pauli suppression of three-body forces in 
dilute neutron matter makes it a good testing ground for chiral effective field theory applied 
to many-nucleon systems. 

The organization of the paper is as follows. We review the lattice interactions contained 
in the leading-order (LO) and next-to-leading-order (NLO) transfer matrices. These transfer 
matrices are rewritten in terms of one-body interactions with auxiliary fields. This allows 
us to simulate the ground state of the many-neutron system using transfer matrix projection 
and hybrid Monte Carlo. The results of the simulations are compared with published results 
for the ground state energy. We also analyze the ground state energy as an expansion near 
the unitarity limit, where the scattering length is infinite and the interactions have negligible 
range. 



II. LATTICE TRANSFER MATRICES WITHOUT AUXILIARY FIELDS 

In [l| we defined the lattice transfer matrix Mloj at leading order and Mnlo2 at next-to- 
leading order. We use the same lattice conventions here and briefly summarize the relevant 
definitions in the appendix. Throughout we use spatial lattice spacing a = (100 MeV)~^ 
and temporal lattice spacing at = (70 MeV)~^. We take for our physical constants m = 
938.92 MeV as the nucleon mass, = 138.08 MeV as the pion mass, = 93 MeV as 
the pion decay constant, and qa = 1-26 as the nucleon axial charge. In ll] we also defined 
lattice actions Mlqi aiid Mnloi • Given the significant computational resources required for 
the Monte Carlo simulations, we use only the Gaussian-smeared actions Mlo2 and Mnlo2- 
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These yield a slightly better description of the 5'-wave interactions expected to be dominant 
in dilute neutron matter. Since we consider only one version of the lattice action here 
we drop the subscript "2" and write Mlq and Mnlo- Tests of model independence using 
different lattice actions and lattice spacings will be pursued in future studies. 
The leading-order lattice transfer matrix is 



Mlo = : exp ^ -i^freett* " ^ /(<? 



+ 



(1) 



Si,S2,I ni,n2 



All of the terms appearing in Eq. ([T]) were defined in |l| and summarized in the appendix. 
The isospin of any two-neutron state is = —1, / = 1. Therefore only the linear combina- 
tion 

C^=i =C + Ci2 (2) 

contributes to systems with only neutrons. As in jl| the coefficient C^^^ is set to —3.414 x 
10-5 MeV-2. 

At next-to-leading order the lattice transfer matrix is 



Ml 



NLO 



exp ^ -iffreettt " ^ ^ /(? 



Cp'^'^^iq)p^'^'^i-q)+Cj.J2pf'^^^pf'^^-^ 



- at [AV + AVi2 + Vg2 + V}2,,2 + ^52^,2 + Vs2j2^g2 + V(,.5)2 + V,2,(,.5)2 + V^-^syk] 



2 ^,2 



+ E E G5.52(-i -n,)pJ;}(n0pSf,(n2) : • 

61,62,/ rai,n2 I 

The NLO corrections to the leading-order contact interactions are 



(3) 



AV = -AC : 2]]p'^^'"(n)p"^''^(n) :, 

n 

Al-,2 = lAC,2:5^pf''^(n)pf''^(n) 



(4) 
(5) 



Again only the 1 = 1 combination contributes to neutron-neutron scattering, 



AC^=^ = AC + AC12. 



(6) 
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There are seven independent NLO contact interactions with two derivatives, 



n,l 



n,I,l 

\i^,q^ = -7;Cs2,p,q2 ■ 2^ Ps,i W^iPs'i W ^ 



n,S,I,l 



n] ■ 



(7) 

(8) 
(9) 
(10) 

(11) 
(12) 



n,I S 



1=1 

{iqxS)-k 



Tlf'^{n)A,pf\n) + Tiffin) A, p^'^\n) 



n,l,S,l' 



n,l,S,l',I 



(13) 



The various static densities, current densities, and symbols A; and S/f, are defined in the 
appendix. The V^^^^"^^).^ interaction is aheady projected onto 1 = 1. The other interactions 
give three independent 1 = 1 coefficients. 



CL ^ = Cg2 + Cl2 q2 



Cs2^g2 = g2 + Cs2j2^g2, 



c, 



1=1 



C, 



(14) 

(15) 
(16) 



'{q-Sr = + Cj2^^g.s)2. 

There are a total of five independent 7 = 1 coefficients at NLO. These were computed 
in [l| using the spherical wall method {3]. The values of the 1=1 coefficients are shown 
in Table [H While the LO terms in the transfer matrix are iterated nonperturbatively, the 
contribution from each NLO interaction is computed using first-order perturbation theory. 
In Fig.[T]the resulting scattering phase shifts for the 7=1 singlet S'-wave and triplet P-waves 
are shown together with partial wave results from . The five arrows show data points used 



TABLE I: Results for 1 = 1 NLO operator coefficients 





-7.7 X 10~'^ MeV-2 


r 


-1.42 X 10-9 MeV"'' 


r'i=i 


-4.53 X 10-1° MeV-"^ 


(-11=1 


-1.80 X 10-10 MeV-'' 


^{iqxS)-k 


9.81 X IQ-ii MeV"^ 



to determine the five 7 = 1 NLO coefficients. There are also four other data points used 
in |l| to determine the four 1 = coefficients, but these are irrelevant for neutron-neutron 
scattering. 

Up until this point our lattice formalism has been constructed for a general system of 
low-energy nucleons. For reasons of numerical efficiency for the Monte Carlo simulation 
we now specialize to the case where all nucleons are neutrons. With this restriction the 
following replacements are possible: 

Cp^'^\q)p^'^\-q) + C,. ^pf ''^(glpf ''^(-g) - C'=' p^' ^\q)p^' '\-q). (17) 
5Z 5Z '^SiS2(^i - ni)pi;}{ni)pi;j{n2) ^ XI 5Z '^5i52(^i - pi;" [ni] p^;" {n2) . 

Si,52,-f ni,n2 5i,52ni,n2 

(18) 

This change has no effect on the interactions between neutrons but leads to the simplified 
transfer matrix, 

Mlo exp | -f7free«* - J] /(g')p'^''"(g1p'^''"(-g1 

+ ll^ E E Gs.S2{n, - n2)pi;\n,)pi;\n2) \ : . (19) 

Sl,S2ni,n2 I 
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FIG. 1: Scattering phase shifts for the 1=1 singlet 5- wave and triplet P- waves versus center- 
of-mass momentum. The five arrows show data points used to determine the five 1=1 NLO 
coefficients. 



At next-to-leading order the simplified neutron transfer matrix is 



NLO 



: exp < -Hirccdt 



2L3 



- at 



C. 



1=1 



C 



V(g.5)2 + y(iqxS)-k 



2 2 

+ J2 Gs,S2{ni-n2)ps;\ni)p''sf{n2) > : . 



(20) 



5i,S2, ni,n2 



In the following we use these simplified forms for the leading-order and next-to-leading-order 
transfer matrices. 
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III. LATTICE TRANSFER MATRICES WITH AUXILIARY FIELDS 



The transfer matrices in Eq. (fT9l) and fl20l) can be rewritten as one-body interactions with 
auxihary fields. This auxihary-field formulation is useful for the many-body simulation. 
The exact equivalence between lattice formalisms with and without auxiliary fields was 
shown in [g], Q]. We summarize the results here. 

In neutron-neutron scattering only the neutral pion contributes to one-pion exchange. Up 
to this point we have been writing the rescaled neutral pion field as vTg. In the following we 
drop the subscript "3" and simply write vr'. Let M^^^\tc', s) be the leading-order auxiliary- 
field transfer matrix at time step rit, 



M'^"'^(7r', s) =: exp <! -Hf^ccOit + ^^^"L y^Ag7r'(n, nt)pp''{n) 



n ) 



(21) 



We can write Mlo as the normalized integral 



[dtt'Ds e-^^"*'-^-*'M("')(7r',s) 

Mlo = ^ r , (22) 

Vn Us e '^'"'^ 



where 5*1"* is the piece of the instantaneous pion action in Eq. flAlSI) containing the neutral 
pion field at time step rit, 

St\7^') = i^7r'(n,n0vr'(n,n0 - ^Y,rr'{n,n,)n\n + l,n,), (23) 

n n,L 

and 5'is''' is the auxiliary- field action at time step nj, 

^ = \Y.s{n,nt)r\n-n')s{n\nt), (24) 

n,n' 

with 

r\n-n') = ^^Y.J^f (25) 
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The NLO interactions require several additional auxiliary fields. Let 

n n,S n,S 

n,S,S' n,l 



n.l.S 



n,l,S 



(26) 



With these extra fields and linear functional U^'^^\e) we define 



M("*)(7r',s,e) =: exp <( -/ffree«, + ^^^A5vr'(n, njp^;'" 



n 

We also define the normalized integral, 



(27) 



(28) 



When all e fields are set to zero we recover Mlc 

M("*)(0) = Mlo. 



(29) 



To first order the NLO interactions in Mnlo can be written as a sum of bilinear derivatives 
of M("*)(£) with respect to the e fields at e = 0, 



Mnlo = Mlo 

2 ^ 

n 



fep(n, n^) 6ep{n, Ut) 



^ 6ep{n,nt) 5ey2p{n,nt) 



]VfK)( 



]VfK)( 



e=0 



+ 



(30) 



e=0 



IV. TRANSFER MATRIX PROJECTION METHOD 



We use the transfer matrix projection method introduced in [8|. First we give a short 
overview using simple continuum notation. Let jvE'^'"'''^^ be a Slater determinant of free- 
particle standing waves in a periodic cube for N neutrons. Let i^LO be the Hamiltonian at 
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leading order, and Hj^lo be the Hamiltonian at next-to-leading order. Let -ffsu(2)y be the 
same as i^LO, but with one-pion exchange turned off by setting qa to zero. As the notation 
suggests, Hsu(2)t/ is invariant under an exact SU(2) intrinsic-spin symmetry. 
Let us define a trial wavefunction 

|vl/(t'))=exp(-/7su(2)/)|^'''^^>. (31) 

In this approach exp (— -f/^su{2)yi') acts as an approximate low-energy filter. In the auxiliary- 
field Monte Carlo calculation this part of the Euclidean time propagation is positive definite 
for any even number of neutrons invariant under the SU(2) intrinsic-spin symmetry 



With this trial wavefunction we define the amplitude, 

Z{t) = exp i-Hi^ot) l^(t')) , (32) 



as well as the transient energy. 



i?Lo(t) = -|[lnZ(t)]. (33) 



In limit of large Euclidean time t we get 

lim ELo(t) = ^o,LO, (34) 



t— >oo 



where -Eq.lo is the energy of the lowest eigenstate |\E'o) of -f^LO with nonzero inner product 
with l^(t'))- 

To compute the expectation value of some operator O we define 

Zo{t) = exp {-H^ot/2) O exp {-H^ot/2) \^{t')) . (35) 

The expectation value of O for |\E'o) is given by the large t limit, 

lim^ = (v&o|0|vPo). (36) 

t^oo Z[t) 

Let -f^NLO be the next-to-leading-order Hamiltonian. Corrections to the energy at next-to- 
leading order can be computed using O = //nlo — -f^LO- Then 

^™ ~^77r ^ -^o.NLO — -E'cLO, (37) 

t^oo Z [t) 

where i?o,NLO is the ground state energy at next-to-leading order. 
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On the lattice we construct I'^it')) using 

|v[/(t')) = (Mgu(2v)^'°|^f^''^>, (38) 

where t' = Lt^at and Lf^ is the number of "outer" time steps. The amphtude Z(t) is 
constructed using 

Z{t) = {^{t')\iM^of^^\^it')), (39) 
where t = L^.ttf and L^. is the number of "inner" time steps. The transient energy 

E^oit + at/2) (40) 

is given by the ratio of the amphtudes for t and t + at, 

(41) 



z{t) • 

The ground state energy -Eq.lo equals the asymptotic limit. 



Eo,Lo= limELo(t + a*/2). (42) 



t— >oo 



For the ground state energy at NLO we compute expectation values of Mnlo and Mlo 
inserted in the middle of a string of LO transfer matrices, 

ZM^^oit) = (^(^')l (Mlo)^*^/' Mnlo {M^of^^/' \^{t')) , (43) 

ZM,oit) = (Mlo)^*'/'Mlo (Mlo)^'-/' \^it')) . (44) 

Clearly Zmlo(^) same as Z{t + at). We use the ratio of amplitudes, 

= 1 - AE^^o{t)at + ■ ■ ■ , (45) 

to define the transient NLO energy correction Ai?NLo(^)- The ellipsis denotes terms which 
are beyond first order in the NLO coefficients. The NLO ground state energy -Eo,nlo is 
calculated using 

-Eo,NLO = -Eo,LO + lim AEo,NLo(i), (46) 

t— >oo 

The Monte Carlo simulation is carried out using the auxiliary-field formulations of the 
transfer matrices. Integrations over auxiliary and pion field configurations are computed 
using hybrid Monte Carlo with endpoint importance sampling. Details of this method can 



be found in the literature 



0,0,3- 
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TABLE II: Monte Carlo results versus exact transfer matrix calculations for the two-neutron spin 
singlet 5 = and spin triplet 5 = 1. 





5 = (MC) 


5 = (exact) 


5 = 1 (MC) 


5=1 (exact) 


^IjO\^ ^ ^tl [ivievj 




— 9 Q9/19 




9^ (T^n 


%IITJ^^ [10^ MeV3] 


4.869(6) 


4.8620 


0.0003(2) 





9(AgNLo(t)) [109 MeV5] 


1.617(3) 


1.6140 


-1.853(4) 


-1.8524 


^(ASnloW) [109 MeV5] 


-4.85(1) 


-4.8419 


-1.851(4) 


-1.8524 


9(AgNLo(t)) [108 MeV^l 


-6.00(1) 


-5.9822 


7.00(2) 


7.0012 


^(AEnloW) [io7 MeV^] 


0.011(7) 





7.8(1) 


7.8743 


A^NLo(i) [MeV] 


-0.0252(3) 


-0.025025 


3.349(7) 


3.3490 



V. PRECISION TESTS 

We use systems of two neutrons to test the auxiliary-field Monte Carlo simulations. We 
calculate the same observables using both the Monte Carlo code and the exact transfer 
matrix without auxiliary fields. We choose a small system so that stochastic errors are 
small enough to expose disagreement at the 0.1% — 1% level. We choose the spatial length 
of lattice to be L = 4 and set the outer time steps Lf^ = 2 and inner time steps Lf. = 4. 

For the first test we choose |\E'^'''^'^^ to be a spin-singlet state built from the Slater deter- 
minant of standing waves and \1IJ2) with 

(0| aij{n) \tpi) oc 5i,o5j,i, (0| aij{n) \tp2) oc (47) 

For the second test we choose a spin-triplet state with standing waves 

(0| aij{n) \^^Jl) oc cos(^)(5,,o5,,i, (0| a,,,(n) 1^2) oc sin(^)<5,,o5,,i. (48) 

Comparisons between Monte Carlo results (MC) and exact transfer matrix calculations 
(exact) are shown in Table [Til The numbers in parentheses are the estimated stochastic 
errors. We see that in each case the agreement is comparable to the estimated stochastic 
error. 
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VI. RESULTS 



We simulate N = 8 and = 12 neutrons on periodic cube lattices with spatial length 
L = 5, 6, 7 lattice units. For each value of and L we set Lt^ = 10 and vary Lt- from 2 to 
12. For we take the Slater determinant formed by standing waves 

(0| ajj(n) \i'2k+i) oc fk{n)8ifi8j^i, (0| aij{n) |^2fc+2) oc fk{n)Si^iSj^i, (49) 

where 

2iTni \ f {r*\ • /'27rni\ 



Mn) = 1, fiin) = cos(^), /2(n) = sin 



/3(n) = cos(^), /4(n) = sin(^), /5(n)=cos(^). (50) 

For = 8 we use k = 0, 1, ■ ■ ■ , 3, and for = 12 we take k = 0, 1, ■ ■ ■ , 5. For each 
value of Lt- a total of about 10® hybrid Monte Carlo trajectories are generated by 1024 
processors, each running completely independent trajectories. Averages and stochastic 
errors are computed by comparing the results of all 1024 processors. 

Let Eq"^*^ be the energy of the ground state for noninteracting neutrons. In Fig. [2] we 
show the dimensionless ratios 

^Lo(t) A^NLo(t) ELo(t) + A^NLo(t) 

Tpfree ' c^free ' rpfrcc ' \ 

versus Euclidean time t for A^ = 8 and L = 5,6, 7. These are labelled using the shorthand 
LO, ANLO, and NLO respectively. The same quantities are shown in Fig. [3] for A^ = 12. 
The lattice calculations for Ai?NLo(^) require an even number of time steps and so fewer 
data points are available. In addition to the Monte Carlo data we plot the asymptotic 
expressions, 

TTirroc 77'iree ' v / 

AENLo(t) ^ -£^0,NLO - -gp.LO ^^-5E t/2 /ggN 
J^frcc ^frec ' ^ ' 

-E'Lo(t) + AENLo(i^) ^ -^0,NLO j^^-^E t , /g^N 

The unknown coefficients A and B, energy gap 6E, and ground state energies -E'o.lo and 
-^^CNLO are determined by least squares fitting. The e~^^'^ dependence in Eq. (15^ comes 
from the contribution of the lowest excited state with energy 6E above the ground state. 
The e"*^^ */^ dependence in Eq. (|53|1 comes from the matrix element of Mnlo between the 

12 



N=S,L = 5 



0.8 



0.6 



« 0.4 



0.2 



1 1 

LO ^ A---; 

ANLO □ i 

NLO I 1 1 

X-/d.f. = LI 



J L_ 



^=8,1 = 6 



0.8 



0.6 



0.4 p 



0.2 



1 1 

LO ; A- -; 

ANLO ■ B i 

NLO I 1 1 

X^/d.f. = LO 



0.1 0.2 0.3 
; (MeV ') 



0.1 0.2 0.3 
f(MeV"') 



N=8,L = 1 



0.8 
„ 0.6 
I 

0.4 
0.2 



1 1 

LO ^ - A : 

ANLO ^ □ ' 

NLO I 1 1 

Z^/d.f. = 1.2 



0.1 0.2 0.3 

t (MeV') 



FIG. 2: Plots of the three energy ratios defined in Eq. (jSip for N = 8 and L = 5,6, 7. These are 
labelled as LO, ANLO, NLO respectively. 
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FIG. 3: Plots of the three energy ratios defined in Eq. (j5ip for = 12 and L = 5, 6, 7. These are 
labelled as LO, ANLO, NLO respectively. 



ground state and the lowest excited state. The reduced chi-square for each fit is shown in 
Fig. [2]and[3l and in each case they are close to 1. 

We calculate the Fermi momentum kp for each neutron spin from the corresponding 
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FIG. 4: Resa 
for FP 1981 
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3, 



APR 1998 



versus Fermi momentum kp. For comparison we show the results 



1$, CMPR v6 and v8' 2003 ^M], SP 2005 Q, and GC 2007 



a. 



density. In our case = Pi = N/ {2L^) and so 



kF = - (Svr^A^) 

Lj 



1/3 



(55) 



In Fig. m we show the results for i?o.NLo/-E'o versus Fermi momentum kp- The error 



bars on i^cNLoZ-^'o represent uncertainties from the asymptotic fits in Eq. ([52^54]). For 
comparison we show other results from the literature: FP 1981 [l^, APR 1998 13], CMPR 
vQ and vS' Q, SP 2005 15|, and GC 2007 [le). We find good agreement near kp = 120 
MeV but there is disagreement whether the slope is positive or negative. 



VII. ANALYSIS AND DISCUSSION 



Neutron matter at fci? ~ 80 MeV is close to the idealized unitarity limit, where the 5- wave 
scattering length is infinite and the range of the interaction is negligible. At lower densities 
corrections due to the scattering length become more important, and at higher densities 
corrections due to the effective range and other effects become important. In the unitarity 
limit the ground state has no dimensionful parameters other than particle density and so 
the ground state energy of the system should obey a simple relation Eq = ^E^'^'^ for some 
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dimensionless constant ^. The universal nature of the unitarity hmit endows it relevance to 
several areas of physics, and in atomic physics the unitarity limit has been studied extensi vely 
with ultracold ^Li and atoms using a magnetic-field Feshbach resonance 



20|. 



Recent experiments for ^ have measured the expansion of ''Li and in the unitarity 



imit released from a harmonic trap. The measured values for ^ are 0.51(4) (2l|, 0.46"*"^^ 



21, and 0.32^}° 



05 



231]. The discrepancy between these measurements and larger values for 



^ reported in earlier experiments 



24 



25 



26| suggests that further work may be needed 



There have been numerous analytic calculations of ^ 



3, 



28 



29 



30 



31 



32 



33 



34 



35 



36|. 



The values for ^ vary roughly from 0.2 to 0.6. Fixed-node Green's function Monte Carlo 



calculations have found ^ to be 0.44(1) j37| and 0.42(1) [38|. An estimate based on Kohn- 
Sham theory for the two-fermion system in a harmonic trap yields a value of 0.42 [39]. There 
have also been simulations of two-component fermions on the lattice in the unitarity limit 



at non-zero temperature. When data are extrapolated to zero temperature the resu 



40| produce a value for ^ similar to the fixed-node results. The same is true for [41 



ts of 



42|, 



though with significant error bars, while calculations by Lee and Schafer 43|, |4J] established 
a bound, 0.07 < ^ < 0.42. 

For finite S-wave scattering length the deviation away from unitarity can be parame- 
terized as 



(56) 



Both f and have been computed using the lattice transfer matrix proiection method 

n n 

discussed above. The results are ^ = 0.25(3) |8| and = 1.0(1) p]. More recent lattice 



calculations find similar values for f an d ^ 



recent literature on the value of 



38 



dij45|, 
48|,M. 



46 



47( 1 . There is general agreement in the 



Further work will be needed to resolve t 



he 



remaining differences between the various calculations for f . We use the values from 

n 

and [6] in our analysis. 

In addition to the corrections at finite scattering length we expect corrections proportional 
to kpTo due to the S'-wave effective range tq. We also expect higher-order corrections away 
from the unitarity limit arising from higher powers of l/lkpao) and kpro, as well as other 
terms associated with the S'-wave shape parameter and triplet P-wave scattering volumes. 
In Fig. [5] we show both Eq^i^q/ Eq'^'^ and Eq^^i^q/Eq^^ versus kp. For comparison we plot 

6 



fikpao) = ^ 



kpQiQ 



(57) 
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FIG. 5: Plot of Eq^i^o/Eq'^^ and Eq^j^i^q / Eq*^^ versus kp- For comparison we plot /{kpao), 
/{kpao) + O.lSfciTTo, and /{kpao) + O.lbkpro + (—1.6 fm^)A:|,. 

with ^ = 0.25, ^1 = 1.0, and neutron scattering length oq = —18.5 fm. From Fig. [5] we see 
that the NLO energy ratio Eq^^i^q / E^'^'^ is approximately described by 

^o,NLo/^o"' - fikpao) + 0. 15fci.ro, (58) 

where ro is the neutron effective range 2.7 fm. The kpro term in Eq. (|58l) can be interpreted 
as the correction due to the neutron effective range. But as noted above there should also be 
corrections from higher powers of l/lkpcio) aiid kprQ and from the S-wave shape parameter 
and triplet P-wave scattering volumes. It is not obvious why these higher-order effects are 
all numerically small at kp ^ as the NLO lattice results suggest. 

In contrast we see deviations beyond the l/{kpaQ) and kpro corrections in the LO lattice 
results. As shown in Fig. [5] the leading-order ratio Eq^i^q / Eq'^^ appears to lie on the curve 

Eo,Lo/^o"' ~ fikpao) + O.lbkpro + (-1.6 fm=^)4. (59) 

We know from the ^Sq phase shifts in Fig. [T]that S-wave scattering for the LO and NLO 
actions are nearly identical. This explains the common coefficient of 0.15 in front of kpro 
for both LO and NLO results. Therefore, the difference between LO and NLO results must 
come from interactions in higher partial waves. 

For the LO action each of triplet P-wave interactions in Fig. [1] are attractive. The 
(—1.6 fm^)kp term in Eq. ( !59|) for the LO action is consistent with the type of correction 



16 



0.4 



3 I 

^Pq removed I \ 1 

^P^ removed --->K---; 

P-, removed B 

. . _ . 



0.3 - 



0.2 - 



3/'— *■ 

all ^ P removed 



.,- (0.18fm^)4 

^ — (0.12 W)4 




0.1 - 



f 




•^■■■(-0.09fm^)4 



-0.2 - 



-0.3 - 



-0.4 - 



"■'■■■(-0.40fm^)4 







50 



100 



150 



200 



(MeV) 



FIG. 6: The change AEq^^i^q / Eq^^ due to removing ^Pq, ^Pi, ^P2, or all triplet P-wave interactions. 

we expect from the negative triplet P-wave scattering volumes. On the other hand, the kp 
correction from P-wave interactions in the NLO action seems to be numerically very small. 
To understand this better we probe the relation between low-energy P-wave interactions 
and the energy ratio Eq^^i^q / Eq"^*^ by varying coefficients of the NLO operators. 

The NLO coefficients in Table [1] were determined by fitting the five data points labelled 
by arrows in Fig. [H We consider four variations of these NLO coefficients. For the first 
variation we set the phase shift for the ^Pq data point to zero while keeping other data points 
the same. For the second variation we zero out the phase shift of the ^Pi data point while 
keeping others the same. For the third variation we zero out only the ^P2 phase shift, and 
for the fourth we zero out all three triplet P- waves. The change Ai^o.NLoZ-^o'^'''' due to each 
of these variations is plotted in Fig. [61 The results show significant cancellation between 
the P-wave contributions. In fact the total contribution from all P-waves is smaller than 
any individual contribution. In Fig. [7] we show Po,nlo / Eq^'^ and the effect of removing all 
triplet P-wave contributions. Both data sets lie approximately on the curve 



We see that the total effect of the triplet P-wave scattering volumes is small due to cancel- 
lations between the J = 0, 1, 2 contributions. 



-Eo.NLo/ E^ 



ifree 




(60) 
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FIG. 7: Plot of Eq^^i^o / Eq versus kp and the effect of removing all triplet P-wave contributions. 



VIII. SUMMARY 

We have discussed lattice simulations of the ground state of dilute neutron matter using 
chiral effective field theory at next-to-leading order. In the first paper the coefficients of the 
next-to-leading-order lattice action were determined by matching nucleon-nucleon scattering 
data for momenta up to the pion mass. In this second paper we used the same lattice action 
to simulate the ground state of up to 12 neutrons in a periodic cube using Monte Carlo for 
the density range from 2% to 8% of normal nuclear density. We found agreement near 
kp = 120 MeV for the ground energy ratio Eq/Eq'^'^ with results in the literature. However 
there is disagreement on whether the ratio is slightly increasing or slightly decreasing with 
kp. 

We analyzed the energy ratio as an expansion about the unitarity limit with corrections 
due to finite scattering length, effective range, and P-wave interactions. We find significant 
cancellation between the various triplet P-wave contributions to the ground state energy. 
We find a good fit to the lattice data using 

Eo,NLo/E'r -^--^^+ 0-15A:Fro, (61) 

with C, = 0.25, ^1 = 1.0. The coefficient in front of kpVQ should be a universal constant and 
therefore measurable in other quantum systems near the unitarity limit. In future studies 
we will consider larger systems of dilute neutron matter and test model independence of 
results in the manner discussed in jl] using several different lattice actions. 
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APPENDIX A: LATTICE ACTION 
1. Notation 

We assume exact isospin symmetry and neglect electromagnetic interactions, n repre- 

— * 

sents integer- valued lattice vectors on a three-dimensional spatial lattice, and p, q, k represent 
integer- valued momentum lattice vectors. / = i, 2, 3 are unit lattice vectors in the spatial 
directions, a is the spatial lattice spacing, and L is the length of the cubic spatial lattice in 
each direction. The lattice time step is at, and rit labels the number of time steps. We 
define at as the ratio between lattice spacings, at — at/ a. Throughout we use dimensionless 
parameters and operators, which correspond with physical values multiplied by the appro- 
priate power of a. Final results are presented in physical units with the corresponding unit 
stated explicitly. 

We use a and to denote annihilation and creation operators. We make explicit all 
spin and isospin indices, 

<20,0 = CtT.P' "^O,! = 0,'\,n, (Al) 

01,0 = ai,p, ai,i = ai,n- (A2) 

The first subscript is for spin and the second subscript is for isospin. We use tj with 
/ = 1,2,3 to represent Pauli matrices acting in isospin space and as with S — 1, 2, 3 to 
represent Pauli matrices acting in spin space. 

We use the eight vertices of a unit cube on the lattice to define spatial derivatives. For 
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each spatial direction / = 1, 2, 3 and any lattice function f{n), let 

^if(^) = l E i-ir^'fin + u), U=v,i + v^2 + vA (A3) 

We also define the double spatial derivative along direction /, 

Vf/(n) = f{n + + f{n - 1) - 2f{n). (A4) 

2. Local densities and currents 



We define the local density, 



which is invariant under Wigner's SU(4) symmetry 501]. Similarly we define the local spin 
density for S = 1,2, 3, 

p1'''iri)= ^lA^)l^s]ii,ai>j{n), (A6) 

i,j,i'=0,l 

isospin density for / = 1,2,3, 

pf '"(^) = Yl Wii'aij'(^), (A7) 

and spin-isospin density for S", J = 1, 2, 3, 

Psfiri)= Y Wjyai',i'(^)- (A8) 

,j'=0,l 

For each static density we also have an associated current density. Similar to the defini- 
tion of the lattice derivative A; in flA3l) . we use the eight vertices of a unit cube, 

z7= z/ii + Z/22 + Z/33, (A9) 

for i^i, 1^2, 1^3 = 0, 1. Let for / = 1, 2, 3 be the result of reflecting the Z*'^- component of 

V about the center of the cube, 

i7(-/) = z7+(l-2z/i)/. (AlO) 

Omitting factors of i and 1/m, we can write the Z^'^-component of the SU(4)-invariant current 
density as 

nf •"(^) = \ Y E ^-^r^"<A^ + ^{-l)y^.{n + y)- (All) 



4 

^1,1^2, 1^3=0,1 1,3=0,1 
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Similarly the Z*^-component of spin current density is 



t^i,i^2,s^3=0,l i,j,i'=0,l 



Z*'^-component of isospin current density is 



and Z*'^-component of spin-isospin current density is 



1^1,1^2,1^3=0,1 i,j,i',j'=0,l 



3. Instantaneous free pion action 



(A12) 



(A13) 



(A14) 



The lattice action for free pions with purely instantaneous propagation is 
S^-rrinj) = at{^ + 3) ^ 7ii{n,nt)ni{n,nt) - at ^ ni{n,nt)7ii{n + i,nt), (A15) 

n,nt,I n,nt,I,l 

where tt/ is the pion field labelled with isospin index /. It is convenient to define a rescaled 
pion field, ttj, 

Trj{n,nt) = y/q^ni{fi,nt), 
q„ = at{ml + 6). 

Then 

'S'ttttW) = ^ E '^i(ri,nt)-Kj(n,nt) - — 7rj(n,nt)nj(n + l,nt). 



n,nt ,1 

In momentum space the action is 



(A16) 
(A17) 

(A18) 



n,nt,I,l 



I,k 



1 _ 

2 Qn 



E-«(^) 



(A19) 



The instantaneous pion correlation function at spatial separation n is 



//^ X //^ x\ J Dtt'j nj{n,ntyj{6,nt) exp [-^^^j 
7rj(n, nt)7rj(0, nt) > = — ^ — -—^ — - — ^ — (no sum on I) 



J Dn'j exp [—5*^ 



(A20) 
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where 



l-!!^E.cos(M.) 



It is useful also to define the two- derivative pion correlator, GsiS2{'^) 



(A21) 



GsiSiin) = (Asyi{n,nt)As2n'i{0,nt)) (no sum on /) 



T6 ^ ^ (-^r^(-^r'''{<i^ + ^-^^^tyAO,n,)). (A22) 



4. Leading-order transfer matrix LO2 



We use the 0(a^)-improved free lattice Hamiltonian, 



n i,j=0,l 



Am 



^ ^ ^ [a^(n)aij(n + /) + 4^.(n)ai,j(n-/) 



n «,j=0,l «=1,2,3 



n i,i=0,l «=1,2,3 



E E E + 30 + «L(^)«i,.(^ - 3/)] . (A23) 



n j,i=0,l /=1,2,3 



The leading-order transfer matrix designated Mloz W\ 



a 



Mlo2 = : exp ^ -i^freettt - ^ X] •^(^ 



+ 



5Z 5Z "^5152(^1 - ^2)Psi;7(^l)P5l'}(^2) 



Si, 82,1 ni,n2 



(A24) 



where the momentum-dependent coefficient function /(g^) is defined as 



/(g^) = /o-^exp 



-b^ (1 - cosg/) 



and the normalization factor /q is determined by the condition 



(A25) 



exp 



-"Ed 



cosg;, 



(A26) 
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The value b — 0.6 gives approximately the correct average effective range for the two S-wSive 
channels when C and C72 are properly tuned. C is the coefficient of the Wigner SU(4)- 
invariant contact interaction and C72 is the coefficient of the isospin-dependent contact 
interaction. For C and Cp we use the values 

C = (3C^=' + C^=°) /4, (A27) 

Cj2 = (C^=i - C^=°) /4, (A28) 
with C^=i = -3.414 X lO'^ MeV'^ and C^=° = -4.780 x lO'^ MeV'^. 
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